home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / c / fftsing.zip / SING.DOC < prev    next >
Text File  |  1992-04-15  |  5KB  |  141 lines
  1. WHY:
  2.     Recently I was asked to Fourier analyze recorded acoustical data.
  3.     I have programs that handle data which does not exceed two 64K
  4.     segments and one program  for extremely long series that will not
  5.     fit in RAM (conventional plus extended). The data I had to analize
  6.     fell between this two extremes and I suspected it could be handled
  7.     easily using huge pointers.
  8.  
  9. DATA SIZE RESTRICTION:
  10.     
  11.     The series has to fit in conventional memory (ie. far heap). Using
  12.     double precision it means somewhere around 70.000 real data points
  13.     (547 K bytes), depending on your system configuration.
  14.  
  15.     If your data exceeds this limit my program TOGAHOCK.PAS might
  16.     be useful or consult the following reference:
  17.  
  18.     Performing Fourier transforms on extremely long data streams
  19.     W. K. Hocking
  20.     Computers in Physics- JAN FEB 1989.
  21.   
  22. TRUNCATION ERRORS ARE A BIG PROBLEM:
  23.     
  24.     First I used two FFT subroutines found in SYMTEL and modified them
  25.     using huge pointers to access the data. Everything worked fine until
  26.     I started transforming series 12K long and above. It took me a while
  27.     to figure out the problem was in the truncation errors when calculating
  28.     the sines and cosines using the library functions.
  29.  
  30. METHOD:
  31.     
  32.     I translated to C, R. C. Singleton's mixed radix fast Fourier transform
  33.     algorithm (see reference in SING.C). His algorithm generates the 
  34.     sines and cosines recursively and corrects for truncation errors.
  35.  
  36. DATA LENGTH
  37.  
  38.     Does not have to be a power of 2 necessarily. The length can contain
  39.     even factors as 2 and 4, and also odd factors as 3,5, 7, 11, etc.
  40.     The algorithm is most efficient if the length is a power of four.
  41.     Data lengths with odd factors of 3 and 5 can be used without a great
  42.     loss    in performance.
  43.  
  44. USAGE TO OBTAIN THE DIRECT TRANSFORM:
  45.         
  46.     In the calling program:
  47.     1-) allocate memory in the far heap for the real and imaginary
  48.         parts (ie. using farcalloc)
  49.     2-) equate two huge pointers to the beginning of the real and
  50.         imaginary parts, respectively
  51.     3-) read data into real and imaginary parts
  52.     4-) call
  53.         sing(huge_points_to_real, huge_points_to_imaginary,
  54.                  order, order, order, -1);
  55.     5-) divide output by order.
  56.  
  57. USAGE TO OBTAIN THE DIRECT TRANSFORM OF REAL DATA USING THE DOUBLING
  58. ALGORITHM: (total_length = 2*half_length)
  59.         
  60.     In the calling program:
  61.     1-) allocate memory in the far heap for the real and imaginary
  62.         parts (ie. using farcalloc), each of length (half_length +1)
  63.     2-) equate two huge pointers to the beginning of the real and
  64.         imaginary parts, respectively
  65.     3-) read real data alternatively into real and imaginary arrays:
  66.           real_array contains  r(0), r(2), r(4), ....
  67.           and imaginary_array  r(1), r(3), r(5).......
  68.         Zero real_array(half_length)and Imaginary_array(half_length)
  69.     4-) call
  70.         sing(huge_points_to_real, huge_points_to_imaginary,
  71.                  half_order,half_order,half_order, -1);
  72.         realtr(huge_points_to_real, huge_points_to_imaginary,
  73.                           half_order, -1);
  74.     5-) divide output by 4*half_order.
  75.  
  76. USAGE TO OBTAIN THE INVERSE TRANSFORM:
  77.         
  78.     In the calling program:
  79.     1-) allocate memory in the far heap for the real and imaginary
  80.         parts (ie. using farcalloc)
  81.     2-) equate two huge pointers to the beginning of the real and
  82.         imaginary parts spectra
  83.     4-) call
  84.         sing(huge_points_to_real, huge_points_to_imaginary,
  85.                  order, order, order, 1);
  86.  
  87. USAGE TO OBTAIN THE INVERSE TRANSFORM WHEN ONLY THE SIMMETRICAL
  88. HALF IS AVAILABLE: (ie, transform of real data using the doubling
  89. algorithm)
  90.         
  91.     In the calling program:
  92.     1-) allocate memory in the far heap for the real and imaginary
  93.         parts (ie. using farcalloc)
  94.     2-) equate two huge pointers to the beginning of the real and
  95.         imaginary parts spectra
  96.     4-) call
  97.         realtr(huge_points_to_real, huge_points_to_imaginary,
  98.                           half_order, +1);
  99.         sing(huge_points_to_real, huge_points_to_imaginary,
  100.                  order, order, order, 1);
  101.     5-) divide output by 4*half_order.            
  102.     6-) the real time series is alternatively in the real and
  103.         imaginary arrays:
  104.           real_array contains  r(0), r(2), r(4), ....
  105.           and imaginary_array  r(1), r(3), r(5).......
  106.  
  107.  
  108. PERFORMANCE :
  109.     Using the doubling algorithm and an AT with 80287 math coprocessor
  110.     
  111.     Total Length       Half Length      Approx time in seconds
  112.  
  113.      2048                 1024                   2
  114.      6144             3072                   8
  115.      8192             4096                  11
  116.     10240             5120                  15
  117.     16384             8192                  24
  118.     20480            10240                  33
  119.     24576            12288                  37
  120.     32768            16384                  48
  121.     40960            20480                  64
  122.     51200            25600                  84
  123.     61440            30720                 112
  124.     65536            32768                 109
  125.     70000            35000                 144
  126.  
  127. VERIFICATION:
  128.  
  129.     Versing.c uses the doubling algorithm. This program generates
  130.     a rectangular pulse, calculates its fourier transform and compares
  131.     it with the known analytic closed form expression. It then
  132.     calculates the inverse transform to get back the original
  133.     rectangular pulse. This is a practical way to asses the
  134.     propagation of truncation errors.
  135.  
  136.  
  137. FINAL COMMENT:
  138.  
  139.     Constructive criticisms and suggestions are welcomed.    I am
  140. willing to adapt this programs to your special needs as long as
  141. I can find free time to do it.